######################################################################################
######################################################################################
#Replication code for:

#Mark T. Buntaine, Brigham Daniels, Colleen Devlin
#Can information outreach increase participation in community-driven development? A field experiment near Bwindi National Park, Uganda
#World Development, 106: 407-421

#Prepared by: Mark Buntaine
#Mark Buntaine contact (as of September 2018): buntaine@bren.ucsb.edu

#Compiled using R Version 3.4.3 (version "Joy in Playing") on Mac running OS X 10.13.6
#+ fig.width=6.5, fig.height=8.25

#Note: all analysis blocks require setup through line 374 run first
######################################################################################
######################################################################################

root <- "~/Google Drive/Bwindi project/Code/WD_Replication"
#Note: all code and data files should be placed in a single directory
#Note: change the file path "root" to the local working directory
setwd(root)

######################################################################################
###Packages
######################################################################################
library(interplot) #Version 0.1.5
library(interflex) #Version 1.0.3
library(lfe) #Version 2.6-2291
library(Lmoments) #Version 1.2-3, for inter.binning.90 function


######################################################################################
###Functions
######################################################################################

mean.imp <- function(data,var,cluster){
  for (i in 1:length(unique(data[,cluster]))){
    imp.value <- mean(data[data[,cluster]==unique(data[,cluster])[i], var], na.rm=T)
    data[data[,cluster]==unique(data[,cluster])[i] & is.na(data[,var]), var] <- imp.value
  }
  return(data)
} #This function automates mean imputation by cluster


felm.ri <- function(formula, dta, treat.var, sims, ...){
  require(lfe)
  
  ate <- coef(felm(formula, data=dta))[1]
  N <- felm(formula, data=dta)$N
  ate.samp.dist <- rep(NA,sims)
  
  for (i in 1:sims){
    dta[,treat.var] <- dta[,ncol(M)+i]
    ate.samp.dist[i] <- coef(felm(formula, data=dta))[1]
  }
  
  p.two.way <- sum(abs(ate)<abs(ate.samp.dist))/sims
  p.one.way.greater <- sum(ate<ate.samp.dist)/sims
  p.one.way.lesser <- sum(ate>ate.samp.dist)/sims
  se <- sd(ate.samp.dist)
  
  coun <- list("ate" = ate, "ate.samp.dist" = ate.samp.dist, "se"=se, "p.two.way" = p.two.way, "p.one.way.greater" = p.one.way.greater, "p.one.way.lesser" = p.one.way.lesser, "N" = N)
  return(coun)
}


lm.ri <- function(formula, dta, treat.var, sims, ...){
  require(lfe)
  
  ate <- coef(lm(formula, data=dta))[2]
  N <- nobs(lm(formula, data=dta))
  ate.samp.dist <- rep(NA,sims)
  
  for (i in 1:sims){
    dta[,treat.var] <- dta[,ncol(M)+i]
    ate.samp.dist[i] <- coef(lm(formula, data=dta))[2]
  }
  
  p.two.way <- sum(abs(ate)<abs(ate.samp.dist))/sims
  p.one.way.greater <- sum(ate<ate.samp.dist)/sims
  p.one.way.lesser <- sum(ate>ate.samp.dist)/sims
  se <- sd(ate.samp.dist)
  
  coun <- list("ate" = ate, "ate.samp.dist" = ate.samp.dist, "se"=se, "p.two.way" = p.two.way, "p.one.way.greater" = p.one.way.greater, "p.one.way.lesser" = p.one.way.lesser, "N" = N)
  return(coun)
}

felm.ri.power <- function(formula, dta, treat.var, sims, ate.bar, outcome.var, te, ...){
  require(lfe)
  
  ate.samp.dist <- rep(NA,sims)
  outcome.store <- dta[,outcome.var]
  
  for (i in 1:sims){
    dta[,treat.var] <- dta[,ncol(M)+i]
    dta[,outcome.var] <- outcome.store
    dta[dta[,treat.var]==1, outcome.var] <- dta[dta[,treat.var]==1, outcome.var] + rbinom(prob = te, size = 1, n=length(dta[dta[,treat.var]==1, outcome.var]))
    dta[,outcome.var] <- ifelse(dta[,outcome.var]>max(outcome.store, na.rm=T), max(outcome.store, na.rm=T), dta[,outcome.var])
    ate.samp.dist[i] <- coef(felm(formula, data=dta))[1]
  }
  
  power <- sum(ate.samp.dist>ate.bar)/sims
  
  coun <- list("ate.samp.dist" = ate.samp.dist, "power"=power)
  return(coun)
}

lm.ri.power <- function(formula, dta, treat.var, sims, ate.bar, outcome.var, te, ...){
  
  ate.samp.dist <- rep(NA,sims)
  outcome.store <- dta[,outcome.var]
  
  for (i in 1:sims){
    dta[,treat.var] <- dta[,ncol(M)+i]
    dta[,outcome.var] <- outcome.store
    dta[dta[,treat.var]==1, outcome.var] <- dta[dta[,treat.var]==1, outcome.var] + rbinom(prob = te, size = 1, n=length(dta[dta[,treat.var]==1, outcome.var]))
    dta[,outcome.var] <- ifelse(dta[,outcome.var]>max(outcome.store, na.rm=T), max(outcome.store, na.rm=T), dta[,outcome.var])
    ate.samp.dist[i] <- coef(lm(formula, data=dta))[2]
  }
  
  power <- sum(ate.samp.dist>ate.bar)/sims
  
  coun <- list("ate.samp.dist" = ate.samp.dist, "power"=power)
  return(coun)
}

source2 <- function(file, start, end, ...) {
  file.lines <- scan(file, what=character(), skip=start-1, nlines=end-start+1, sep='\n')
  file.lines.collapsed <- paste(file.lines, collapse='\n')
  source(textConnection(file.lines.collapsed), ...)
} #Function to source parts of main replication file

source('inter.binning.90.R')
source('multiplot.R')

#Function to reset par() show that blocks can be run sequentially and without a plotting device
default.par <- par()


######################################################################################
###Data input
######################################################################################

B <- read.csv("baseline.csv", stringsAsFactors=FALSE)
E <- read.csv("endline.csv", stringsAsFactors=FALSE)
endline.subjects <- read.csv("endline_call_list.csv", row.names=1, stringsAsFactors=FALSE)[,1:4]
recruited.subjects <- read.csv("recruited_subjects.csv", stringsAsFactors=FALSE)
smsid <- read.csv("sms_received_list.csv", stringsAsFactors=FALSE)
Responsiveness <- read.csv("digital_responses.csv", stringsAsFactors=FALSE)
con.mat <- read.csv("contiguity_mat1.csv", stringsAsFactors=FALSE)
con2.mat <- read.csv("contiguity_mat2.csv", stringsAsFactors=FALSE)
con3.mat <- read.csv("contiguity_mat3.csv", stringsAsFactors=FALSE)
#Note the "con.mat" files have village "Kiziiba" manually removed because merged with "Buguro" at cleaning phase
spill.dist <- read.csv("spillover_matrix.csv", row.names=1)
loc.master <- read.csv("location_master.csv", stringsAsFactors=FALSE)
lc1.planning <- read.csv("planning_meeting_attendance.csv", stringsAsFactors=FALSE)
sbj_RI <- read.csv("assignment_ri.csv", stringsAsFactors=FALSE)



######################################################################################
###Data setup
######################################################################################

#Removing duplicate subject.id's in endline.subjects
endline.dup <- endline.subjects$Subject_ID_number[duplicated(endline.subjects$Subject_ID_number)]
see <- endline.subjects[endline.subjects$Subject_ID_number %in% endline.dup,]
see <- see[order(see$Subject_ID_number),] #Visual inspection: all subject id's are exact duplicates
M <- endline.subjects[!duplicated(endline.subjects$Subject_ID_number),] #Removing duplicated subject id's
length(unique(M$Subject_ID_number)) #This is the merge point, as all subject id's are unique

#Preparing baseline merge
B <- subset(B, B_subject_refuse_to_participate=="no") #Removing subjects who refused, no data
length(unique(B$Subject_ID_number)) #There are several duplicate subject id's in the raw baseline file
B.duplicated.list <- B$Subject_ID_number[duplicated(B$Subject_ID_number)] #List of duplicate subject id's in raw endline file
see <- B[B$Subject_ID_number %in% B.duplicated.list,]
see <- see[order(see$Subject_ID_number),]
B <- B[!(B$Subject_ID_number %in% B.duplicated.list),]
M <- merge(M,B, by="Subject_ID_number", all.x=TRUE) #still 1723 observations

#Preparing endline merge
E <- subset(E, E_1_Did_the_subject_refuse_to_p=="no") #Removing subjects who refused, no data
E.duplicated.list <- E$Subject_ID_number[duplicated(E$Subject_ID_number)] #List of duplicate subject id's in raw endline file
see <- E[E$Subject_ID_number %in% E.duplicated.list,]
see <- see[order(see$Subject_ID_number),]
E <- E[!(E$Subject_ID_number %in% E.duplicated.list),]
M <- merge(M,E, by="Subject_ID_number", all.x=TRUE) #still 1723 observations

#Merging treatment status and block from the assignment file
treat.dta <- recruited.subjects[,c(1,10,11,12)] #Getting subcounty, treatment status, and block
M <- merge(M,treat.dta, by="Subject_ID_number", all.x=TRUE) #still 1723 observations
M$subcounty.blocking <- ifelse(M$subcounty.blocking=="Rutenga", "Singles", M$subcounty.blocking) #Because of exclusion, Rutenga subcounty has only one village, moving to "Singles" block

#TextIt data
m7<-table(smsid$Subject_ID_number)
m7<-data.frame(m7)
names(m7)<-c("Subject_ID_number","participate.m7")
M <- merge(M,m7, by="Subject_ID_number", all.x=TRUE) #still 1723 observations

#Location cleaning
villages <- loc.master$village.updated
for (i in 1:length(villages)){
  villages[i] <- ifelse(loc.master$village.cleaned[i]=="", loc.master$village.updated[i], loc.master$village.cleaned[i])
}
villages[duplicated(villages)] #Checking for errors; looks good though some original errors resulted in villages with mixed treatment assignment (Kajumo, Omumbuga)

see <- subset(loc.master,excluded==0)
for (i in 1:length(villages)){
  villages[i] <- ifelse(see$village.cleaned[i]=="", see$village.updated[i], see$village.cleaned[i])
}
villages[duplicated(villages)] #Checking for errors; looks good though some original errors resulted in villages with mixed treatment assignment (Kajumo, Omumbuga)
see$final.id <- ifelse(see$location.id.cleaned=="",see$location.id,see$location.id.cleaned)
unique(see$final.id)

#unique(M$Updated.Village)[!(unique(M$Updated.Village) %in% villages)]
M$Updated.Village <- ifelse(M$Updated.Village=="Katooma","Katoma",M$Updated.Village)
M$Updated.Village <- ifelse(M$Updated.Village=="Kijumwa","Kijuma",M$Updated.Village)
M$Updated.Village <- ifelse(M$Updated.Village=="Kyogo (Kabale)","Kyogo (Rubanda)",M$Updated.Village)
M$Updated.Village <- ifelse(M$Updated.Village=="Mushija","Musheija",M$Updated.Village)
M$Updated.Village <- ifelse(M$Updated.Village=="Ndeego","Ndego",M$Updated.Village)

#Adding variables to probe spillover
spill.village.list <- names(spill.dist)
spill.village.list[!(spill.village.list %in% recruited.subjects$updated.village)] #list in spillover matrix that don't match

#con.village.list <- con.mat$Village..m.
#con.village.list[!(con.village.list %in% M$Updated.Village)] #list in contiguity matrix that don't match M village names
#con.village.list[!(con.village.list %in% loc.master$village.updated)] #list in contiguity matrix that don't match location master original names

#Cleaning the con.mat file for village names, see email thread "village names in matrices" from Jan-2017
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kashija/ Rubuguri","Kashija",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Higabiro/ Rubuguri","Higabiro (Rubuguri)",con.mat$Village..m.)
#Higabiro (Rubuguri): this is an RS village, but we do not have any subjects from this village
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Nteeko","Nteko",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Nyakabungo( Rubimbwa)","Nyakabungo (Rubimbwa)",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Nyakabungo(Bushura)","Nyakabungo (Bushura)",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Ishaya (Mucusye)","Ishaya",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kiziba","Kiziiba",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Bushegyenyi","Bushegenyi",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kyibingo","Kibingo",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kitahulira","Kitahurira (Hakikome)",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Nyakabungo","Nyakabungo (Kashekyera)",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Buzaniza","Buzaniro",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Mburamaizi","Mburameizi",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Bishanyu","Bishayu",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Katooma","Katoma",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Rwensaniro","Rwensanziro",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Rugando","Rugandu",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Ntunagamo","Ntungamo",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Nyakaranaga","Nyakaranga",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kacwamuhoro","Kachwamuhoro",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kitahurira","Kitahurira (Kashasha)",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Ndeego","Ndego",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="kagogo","Kagogo",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Byamihanda","Ryamihanda",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Mushija","Musheija",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kijumwa","Kijuma",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kiziiba","Bugoro",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kyogo (Kabale)","Kyogo (Rubanda)",con.mat$Village..m.)

#Cleaning the con2.mat file for village names, see email thread "village names in matrices" from Jan-2017
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kashija/ Rubuguri","Kashija",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Higabiro/ Rubuguri","Higabiro (Rubuguri)",con2.mat$Village..m.)
#Higabiro (Rubuguri): this is an RS village, but we do not have any subjects from this village
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Nteeko","Nteko",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Nyakabungo( Rubimbwa)","Nyakabungo (Rubimbwa)",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Nyakabungo(Bushura)","Nyakabungo (Bushura)",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Ishaya (Mucusye)","Ishaya",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kiziba","Kiziiba",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Bushegyenyi","Bushegenyi",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kyibingo","Kibingo",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kitahulira","Kitahurira (Hakikome)",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Nyakabungo","Nyakabungo (Kashekyera)",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Buzaniza","Buzaniro",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Mburamaizi","Mburameizi",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Bishanyu","Bishayu",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Katooma","Katoma",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Rwensaniro","Rwensanziro",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Rugando","Rugandu",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Ntunagamo","Ntungamo",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Nyakaranaga","Nyakaranga",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kacwamuhoro","Kachwamuhoro",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kitahurira","Kitahurira (Kashasha)",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Ndeego","Ndego",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="kagogo","Kagogo",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Byamihanda","Ryamihanda",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Mushija","Musheija",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kijumwa","Kijuma",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kiziiba","Bugoro",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kyogo (Kabale)","Kyogo (Rubanda)",con2.mat$Village..m.)

#Cleaning the con3.mat file for village names, see email thread "village names in matrices" from Jan-2017
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kashija/ Rubuguri","Kashija",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Higabiro/ Rubuguri","Higabiro (Rubuguri)",con3.mat$Village..m.)
#Higabiro (Rubuguri): this is an RS village, but we do not have any subjects from this village
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Nteeko","Nteko",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Nyakabungo( Rubimbwa)","Nyakabungo (Rubimbwa)",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Nyakabungo(Bushura)","Nyakabungo (Bushura)",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Ishaya (Mucusye)","Ishaya",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kiziba","Kiziiba",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Bushegyenyi","Bushegenyi",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kyibingo","Kibingo",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kitahulira","Kitahurira (Hakikome)",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Nyakabungo","Nyakabungo (Kashekyera)",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Buzaniza","Buzaniro",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Mburamaizi","Mburameizi",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Bishanyu","Bishayu",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Katooma","Katoma",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Rwensaniro","Rwensanziro",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Rugando","Rugandu",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Ntunagamo","Ntungamo",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Nyakaranaga","Nyakaranga",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kacwamuhoro","Kachwamuhoro",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kitahurira","Kitahurira (Kashasha)",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Ndeego","Ndego",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="kagogo","Kagogo",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Byamihanda","Ryamihanda",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Mushija","Musheija",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kijumwa","Kijuma",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kiziiba","Bugoro",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kyogo (Kabale)","Kyogo (Rubanda)",con3.mat$Village..m.)

#Cleaning loc.master names into one vector
loc.master$village.final <- ifelse(loc.master$village.cleaned=="", loc.master$village.updated, loc.master$village.cleaned)
loc.master$village.final <- ifelse(loc.master$village.final=="Kikangaga ", "Kikangaga", loc.master$village.final)

#Checking the location.id cleaning in all files
con.mat$Village..m.[!(con.mat$Village..m. %in% loc.master$village.final)] #Notes: only Higabiro (Rubuguri), no subjects in study
loc.master$village.final[!(loc.master$village.final %in% con.mat$Village..m.)] #Notes: "Kinaba" "Kinaba" "Kyogo"  "Rukere" all excluded
unique(M$Updated.Village) %in% loc.master$village.final #Note: all appear in cleaned loc.master
unique(M$Updated.Village)[!(unique(M$Updated.Village) %in% con.mat$Village..m.)] #"Kinaba" does not specify which one, remove Kinaba from dataframe?

#Correcting Kinaba village names (redistricted), for those we were able to reach, see "Corrected villages.xls" and email thread "Bwindi Village Confusion"
M$Updated.Village[M$Subject_ID_number %in% c(2581,2583,2587)] <- "Nyabiha"
M$Updated.Village[M$Subject_ID_number %in% c(2518)] <- "Kyabworo"
M <- M[M$Updated.Village!="Kinaba",] #Removing "Kinaba" from dataframe per exclusion from revenue sharing during '16/'17 round

#Creating spillover indicator in M
treated.villages <- loc.master$village.final[loc.master$treat==1]
#Note: rows/columns in con.mat must be ordered identically, code below checks
#col.n <- colnames(con.mat)[-1]
#ex <- dataframe(con.mat$Village..m.,col.n)
#con.mat$Village..m.[48] <- "Kiziba" #Correction discussed in 2/20/17 email to Colleen

for (i in 1:nrow(M)){
  row.test <- con.mat[M$Updated.Village[i]==con.mat$Village..m.,-1]==1
  row.test2 <- con2.mat[M$Updated.Village[i]==con2.mat$Village..m.,-1]==1
  row.test3 <- con3.mat[M$Updated.Village[i]==con3.mat$Village..m.,-1]==1
  con.villages <- con.mat$Village..m.[row.test]
  con2.villages <- con2.mat$Village..m.[row.test2]
  con3.villages <- con3.mat$Village..m.[row.test3]
  M$S1[i] <- sum(con.villages %in% treated.villages)
  M$S2[i] <- sum(con2.villages %in% treated.villages)
  M$S3[i] <- sum(con3.villages %in% treated.villages)
  
}
M$S1 <- as.factor(M$S1)
M$S2 <- as.factor(M$S2)
M$S3 <- as.factor(M$S3)

#Treatment Density (number of subjects per village)
for (i in 1:nrow(M)){
  sub <- subset(M, Updated.Village==M$Updated.Village[i])
  M$subjects.per.village[i] <- nrow(sub)
}

#Attrition
M$attr <- ifelse(is.na(M$E_id), 1, 0)

#Setting up RI object
source('outcome_var_mean_imp_setup.R')
sbj_RI <- sbj_RI[,c(1,5:ncol(sbj_RI))] #selecting only the treatment status to avoid introducing duplicate data to M
M.ri <- merge(M, sbj_RI, by="Subject_ID_number", all.x=TRUE)

######################################################################################
###Tables D1-3. Results for F-Test of Improved Model Fit with Spillover Indicators
#Note: models presented in this section exclude most covariates
######################################################################################
#Note: tables formatted manually based on values displayed below

source('outcome_var_mean_imp_setup.R')

###H7.1: Satisfaction with RS Info from UWA
mod.h7 <- lm(e.knowledge.b7 ~ treat + b.knowledge.b7 + subcounty.blocking, data=M)

mod.h7.spill <- lm(e.knowledge.b7 ~ treat + S1 + b.knowledge.b7 + subcounty.blocking, data=M)

mod.h7.spill2 <- lm(e.knowledge.b7 ~ treat + S2 + b.knowledge.b7 + subcounty.blocking, data=M)

mod.h7.spill3 <- lm(e.knowledge.b7 ~ treat + S3 + b.knowledge.b7 + subcounty.blocking, data=M)

###H7.2: Know How RS Works Generally
mod.h7.2 <- lm(e.knowledge.b8 ~ treat + b.knowledge.b8 + subcounty.blocking, data=M)

mod.h7.2.spill <- lm(e.knowledge.b8 ~ treat + S1 + b.knowledge.b8 + subcounty.blocking, data=M)

mod.h7.2.spill2 <- lm(e.knowledge.b8 ~ treat + S2 + b.knowledge.b8 + subcounty.blocking, data=M)

mod.h7.2.spill3 <- lm(e.knowledge.b8 ~ treat + S3 + b.knowledge.b8 + subcounty.blocking, data=M)

###H7.3: How RS Works in my Village
mod.h7.3 <- lm(e.knowledge.b9 ~ treat + b.knowledge.b9 + subcounty.blocking, data=M)

mod.h7.3.spill <- lm(e.knowledge.b9 ~ treat + S1 + b.knowledge.b9 + subcounty.blocking, data=M)

mod.h7.3.spill2 <- lm(e.knowledge.b9 ~ treat + S2 + b.knowledge.b9 + subcounty.blocking, data=M)

mod.h7.3.spill3 <- lm(e.knowledge.b9 ~ treat + S3 + b.knowledge.b9 + subcounty.blocking, data=M)

###H8.1: Opportunities to Participate in Planning
mod.h8.1 <- lm(e.participate.b10 ~ treat + b.participate.b10 + subcounty.blocking, data=M)

mod.h8.1.spill <- lm(e.participate.b10 ~ treat + S1 + b.participate.b10 + subcounty.blocking, data=M)

mod.h8.1.spill2 <- lm(e.participate.b10 ~ treat + S2 + b.participate.b10 + subcounty.blocking, data=M)

mod.h8.1.spill3 <- lm(e.participate.b10 ~ treat + S3 + b.participate.b10 + subcounty.blocking, data=M)

###H8.3: Participate in RS meetings in past months
mod.h8.3 <- lm(e.participate.e4 ~ treat + subcounty.blocking, data=M)

mod.h8.3.spill <- lm(e.participate.e4 ~ treat + S1 + subcounty.blocking, data=M)

mod.h8.3.spill2 <- lm(e.participate.e4 ~ treat + S2 + subcounty.blocking, data=M)

mod.h8.3.spill3 <- lm(e.participate.e4 ~ treat + S3 + subcounty.blocking, data=M)

###H9.1: Satisfaction with opportunities to communicate with UWA about RS
mod.h9.1 <- lm(e.participate.b11 ~ treat + b.participate.b11 + subcounty.blocking, data=M)

mod.h9.1.spill <- lm(e.participate.b11 ~ treat + S1 + b.participate.b11 + subcounty.blocking, data=M)

mod.h9.1.spill2 <- lm(e.participate.b11 ~ treat + S2 + b.participate.b11 + subcounty.blocking, data=M)

mod.h9.1.spill3 <- lm(e.participate.b11 ~ treat + S3 + b.participate.b11 + subcounty.blocking, data=M)

###H9.2: Know Right Person to Contact
mod.h9.2 <- lm(e.participate.b12 ~ treat + b.participate.b12 + subcounty.blocking, data=M)

mod.h9.2.spill <- lm(e.participate.b12 ~ treat + S1 + b.participate.b12 + subcounty.blocking, data=M)

mod.h9.2.spill2 <- lm(e.participate.b12 ~ treat + S2 + b.participate.b12 + subcounty.blocking, data=M)

mod.h9.2.spill3 <- lm(e.participate.b12 ~ treat + S3 + b.participate.b12 + subcounty.blocking, data=M)

###H10: Satisfaction with Park Management
nod.h10mgmt <- lm(e.satisfaction.b2 ~ treat + b.satisfaction.b2 + subcounty.blocking, data=M)

nod.h10mgmt.spill <- lm(e.satisfaction.b2 ~ treat + S1 + b.satisfaction.b2 + subcounty.blocking, data=M)

nod.h10mgmt.spill2 <- lm(e.satisfaction.b2 ~ treat + S2 + b.satisfaction.b2 + subcounty.blocking, data=M)

nod.h10mgmt.spill3 <- lm(e.satisfaction.b2 ~ treat + S3 + b.satisfaction.b2 + subcounty.blocking, data=M)

###H11: Satisfaction with RS
nod.h11 <- lm(e.satisfaction.b3 ~ treat + b.satisfaction.b3 + subcounty.blocking, data=M)

nod.h11.spill <- lm(e.satisfaction.b3 ~ treat + S1 + b.satisfaction.b3 + subcounty.blocking, data=M)

nod.h11.spill2 <- lm(e.satisfaction.b3 ~ treat + S2 + b.satisfaction.b3 + subcounty.blocking, data=M)

nod.h11.spill3 <- lm(e.satisfaction.b3 ~ treat + S3 + b.satisfaction.b3 + subcounty.blocking, data=M)

###H13: Support for Conservation (Importance of Protecting Bwindi)
nod.h13 <- lm(e.conservation.b6 ~ treat + b.conservation.b6 + subcounty.blocking, data=M)

nod.h13.spill <- lm(e.conservation.b6 ~ treat + S1 + b.conservation.b6 + subcounty.blocking, data=M)

nod.h13.spill2 <- lm(e.conservation.b6 ~ treat + S2 + b.conservation.b6 + subcounty.blocking, data=M)

nod.h13.spill3 <- lm(e.conservation.b6 ~ treat + S3 + b.conservation.b6 + subcounty.blocking, data=M)

#Table D1
anova(nod.h13, nod.h13.spill)
anova(nod.h11, nod.h11.spill)
anova(nod.h10mgmt, nod.h10mgmt.spill)
anova(mod.h9.2, mod.h9.2.spill)
anova(mod.h9.1, mod.h9.1.spill)
anova(mod.h8.3, mod.h8.3.spill)
anova(mod.h8.1, mod.h8.1.spill)
anova(mod.h7.3, mod.h7.3.spill)
anova(mod.h7.2, mod.h7.2.spill)
anova(mod.h7, mod.h7.spill)

#Table D2
anova(nod.h13, nod.h13.spill2)
anova(nod.h11, nod.h11.spill2)
anova(nod.h10mgmt, nod.h10mgmt.spill2)
anova(mod.h9.2, mod.h9.2.spill2)
anova(mod.h9.1, mod.h9.1.spill2)
anova(mod.h8.3, mod.h8.3.spill2)
anova(mod.h8.1, mod.h8.1.spill2)
anova(mod.h7.3, mod.h7.3.spill2)
anova(mod.h7.2, mod.h7.2.spill2)
anova(mod.h7, mod.h7.spill2)

#Table D3
anova(nod.h13, nod.h13.spill3)
anova(nod.h11, nod.h11.spill3)
anova(nod.h10mgmt, nod.h10mgmt.spill3)
anova(mod.h9.2, mod.h9.2.spill3)
anova(mod.h9.1, mod.h9.1.spill3)
anova(mod.h8.3, mod.h8.3.spill3)
anova(mod.h8.1, mod.h8.1.spill3)
anova(mod.h7.3, mod.h7.3.spill3)
anova(mod.h7.2, mod.h7.2.spill3)
anova(mod.h7, mod.h7.spill3)


######################################################################################
###Figures D1-3. Results for F-Test of Improved Model Fit with Spillover Indicators
#Note: models presented in this section exclude most covariates
######################################################################################

#Figure D1: modifying only models where 2-village spillover indicator improves model fit (custom style)
names <- c("Satisfied with Information from UWA about RS","Know How RS Works Generally","Know How RS Works in my Village","Opportunities to Participate in RS","Participated in RS meeting (past month)","Opportunities to Communicate with UWA about RS","Know Right Person to Contact about RS","Satisfied with Park Management","Satisfied with Revenue Sharing","Importance of Protecting Bwindi")
ates <- c(mod.h7$coef[2],mod.h7.2$coef[2],mod.h7.3$coef[2],mod.h8.1$coef[2],mod.h8.3.spill$coef[2],mod.h9.1$coef[2],mod.h9.2$coef[2],nod.h10mgmt$coef[2],nod.h11.spill$coef[2],nod.h13$coef[2])
ses <- c(coef(summary(mod.h7))[2,2],coef(summary(mod.h7.2))[2,2],coef(summary(mod.h7.3))[2,2],coef(summary(mod.h8.1))[2,2],coef(summary(mod.h8.3.spill))[2,2],coef(summary(mod.h9.1))[2,2],coef(summary(mod.h9.2))[2,2],coef(summary(nod.h10mgmt))[2,2],coef(summary(nod.h11.spill))[2,2],coef(summary(nod.h13))[2,2])
dta.coef <- data.frame(names,ates,ses)
dta.coef$ylo <- dta.coef$ates - 1.645*dta.coef$ses
dta.coef$yhi <- dta.coef$ates + 1.645*dta.coef$ses
dta.coef$names <- factor(dta.coef$names, levels=dta.coef$names)

spill = ggplot(dta.coef, aes(x=names, y=ates, ymin=ylo, ymax=yhi)) + geom_errorbar(colour="blue", width=0.2, size=1.4) + 
  theme_bw(base_size=16) + coord_flip() + geom_hline(aes(yintercept=0), lty=2) + ylab("Treatment Effects on Survey Items") + xlab("") + geom_point(size=4, colour="blue") +
  theme(axis.title = element_text(size=14))
spill 


#Figure D2: modifying only models where 4-village spillover indicator improves model fit (custom style)
names <- c("Satisfied with Information from UWA about RS","Know How RS Works Generally","Know How RS Works in my Village","Opportunities to Participate in RS","Participated in RS meeting (past month)","Opportunities to Communicate with UWA about RS","Know Right Person to Contact about RS","Satisfied with Park Management","Satisfied with Revenue Sharing","Importance of Protecting Bwindi")
ates <- c(mod.h7$coef[2],mod.h7.2$coef[2],mod.h7.3$coef[2],mod.h8.1$coef[2],mod.h8.3$coef[2],mod.h9.1$coef[2],mod.h9.2$coef[2],nod.h10mgmt$coef[2],nod.h11.spill2$coef[2],nod.h13$coef[2])
ses <- c(coef(summary(mod.h7))[2,2],coef(summary(mod.h7.2))[2,2],coef(summary(mod.h7.3))[2,2],coef(summary(mod.h8.1))[2,2],coef(summary(mod.h8.3))[2,2],coef(summary(mod.h9.1))[2,2],coef(summary(mod.h9.2))[2,2],coef(summary(nod.h10mgmt))[2,2],coef(summary(nod.h11.spill2))[2,2],coef(summary(nod.h13))[2,2])
dta.coef <- data.frame(names,ates,ses)
dta.coef$ylo <- dta.coef$ates - 1.645*dta.coef$ses
dta.coef$yhi <- dta.coef$ates + 1.645*dta.coef$ses
dta.coef$names <- factor(dta.coef$names, levels=dta.coef$names)

spill2 = ggplot(dta.coef, aes(x=names, y=ates, ymin=ylo, ymax=yhi)) + geom_errorbar(colour="blue", width=0.2, size=1.4) + 
  theme_bw(base_size=16) + coord_flip() + geom_hline(aes(yintercept=0), lty=2) + ylab("Treatment Effects on Survey Items") + xlab("") + geom_point(size=4, colour="blue") +
  theme(axis.title = element_text(size=14))
spill2 #Copied at 800x420


#Figure D3: modifying only models where 6-village spillover indicator improves model fit (custom style)
names <- c("Satisfied with Information from UWA about RS","Know How RS Works Generally","Know How RS Works in my Village","Opportunities to Participate in RS","Participated in RS meeting (past month)","Opportunities to Communicate with UWA about RS","Know Right Person to Contact about RS","Satisfied with Park Management","Satisfied with Revenue Sharing","Importance of Protecting Bwindi")
ates <- c(mod.h7$coef[2],mod.h7.2$coef[2],mod.h7.3$coef[2],mod.h8.1$coef[2],mod.h8.3$coef[2],mod.h9.1.spill3$coef[2],mod.h9.2$coef[2],nod.h10mgmt$coef[2],nod.h11.spill3$coef[2],nod.h13$coef[2])
ses <- c(coef(summary(mod.h7))[2,2],coef(summary(mod.h7.2))[2,2],coef(summary(mod.h7.3))[2,2],coef(summary(mod.h8.1))[2,2],coef(summary(mod.h8.3))[2,2],coef(summary(mod.h9.1.spill3))[2,2],coef(summary(mod.h9.2))[2,2],coef(summary(nod.h10mgmt))[2,2],coef(summary(nod.h11.spill3))[2,2],coef(summary(nod.h13))[2,2])
dta.coef <- data.frame(names,ates,ses)
dta.coef$ylo <- dta.coef$ates - 1.645*dta.coef$ses
dta.coef$yhi <- dta.coef$ates + 1.645*dta.coef$ses
dta.coef$names <- factor(dta.coef$names, levels=dta.coef$names)

spill3 = ggplot(dta.coef, aes(x=names, y=ates, ymin=ylo, ymax=yhi)) + geom_errorbar(colour="blue", width=0.2, size=1.4) + 
  theme_bw(base_size=16) + coord_flip() + geom_hline(aes(yintercept=0), lty=2) + ylab("Treatment Effects on Survey Items") + xlab("") + geom_point(size=4, colour="blue") +
  theme(axis.title = element_text(size=14))
spill3 #Copied at 800x420


######################################################################################
###Table D4. Treatment effects under the base specification and adjusted for spillover
######################################################################################
#Note: this block produces the values for table; formatting done separately
source('ols.R') #Function to compute OLS w/ clustered SEs

###H7.2: Know How RS Works Generally
mod.h7.2 <- ols(e.knowledge.b8 ~ treat + b.knowledge.b8 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster = "subcounty.blocking")
mod.h7.2

mod.h7.2.spill <- ols(e.knowledge.b8 ~ treat + S1 + b.knowledge.b8  + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster = "subcounty.blocking")
mod.h7.2.spill

mod.h7.2.spill2 <- ols(e.knowledge.b8 ~ treat + S2 + b.knowledge.b8  + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster = "subcounty.blocking")
mod.h7.2.spill2

mod.h7.2.spill3 <- ols(e.knowledge.b8 ~ treat + S3 + b.knowledge.b8  + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster = "subcounty.blocking")
mod.h7.2.spill3

###H7.3: How RS Works in my Village
mod.h7.3 <- ols(e.knowledge.b9 ~ treat + b.knowledge.b9 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h7.3

mod.h7.3.spill <- ols(e.knowledge.b9 ~ treat + S1 + b.knowledge.b9  + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h7.3.spill

mod.h7.3.spill2 <- ols(e.knowledge.b9 ~ treat + S2 + b.knowledge.b9  + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h7.3.spill2

mod.h7.3.spill3 <- ols(e.knowledge.b9 ~ treat + S3 + b.knowledge.b9  + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h7.3.spill3

###H9.2: Know Right Person to Contact
mod.h9.2 <- ols(e.participate.b12 ~ treat + b.participate.b12 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h9.2
#To Do: the SI entries in this cell are incorrect, fix

mod.h9.2.spill <- ols(e.participate.b12 ~ treat + S1 + b.participate.b12 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster = "subcounty.blocking")
mod.h9.2.spill

mod.h9.2.spill2 <- ols(e.participate.b12 ~ treat + S2 + b.participate.b12 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster = "subcounty.blocking")
mod.h9.2.spill2

mod.h9.2.spill3 <- ols(e.participate.b12 ~ treat + S3 + b.participate.b12 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster = "subcounty.blocking")
mod.h9.2.spill3

###H9.1: Satisfaction with opportunities to communicate with UWA about RS
mod.h9.1 <- ols(e.participate.b11 ~ treat + b.participate.b11 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h9.1

mod.h9.1.spill <- ols(e.participate.b11 ~ treat + S1 + b.participate.b11 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster = "subcounty.blocking")
mod.h9.1.spill

mod.h9.1.spill2 <- ols(e.participate.b11 ~ treat + S2 + b.participate.b11 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster = "subcounty.blocking")
mod.h9.1.spill2

mod.h9.1.spill3 <- ols(e.participate.b11 ~ treat + S3 + b.participate.b11 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster = "subcounty.blocking")
mod.h9.1.spill3

###H8.1: Opportunities to Participate in Planning
mod.h8.1 <- ols(e.participate.b10 ~ treat + b.participate.b10 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h8.1

mod.h8.1.spill <- ols(e.participate.b10 ~ treat + S1 + b.participate.b10 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h8.1.spill

mod.h8.1.spill2 <- ols(e.participate.b10 ~ treat + S2 + b.participate.b10 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h8.1.spill2

mod.h8.1.spill3 <- ols(e.participate.b10 ~ treat + S3 + b.participate.b10 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h8.1.spill3

###H8.3: Participate in RS meetings in past months
mod.h8.3 <- ols(e.participate.e4 ~ treat + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h8.3

mod.h8.3.spill <- ols(e.participate.e4 ~ treat + S1 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h8.3.spill

mod.h8.3.spill2 <- ols(e.participate.e4 ~ treat + S2 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h8.3.spill2

mod.h8.3.spill3 <- ols(e.participate.e4 ~ treat + S3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h8.3.spill3

###H7.1: Satisfaction with RS Info from UWA
#Note: this row is labeled incorrectly in the online SI
mod.h7 <- ols(e.knowledge.b7 ~ treat + b.knowledge.b7 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h7

mod.h7.spill <- ols(e.knowledge.b7 ~ treat + S1 + b.knowledge.b7 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h7.spill

mod.h7.spill2 <- ols(e.knowledge.b7 ~ treat + S2 + b.knowledge.b7 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h7.spill2

mod.h7.spill3 <- ols(e.knowledge.b7 ~ treat + S3 + b.knowledge.b7 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subcounty.blocking, data=M, cluster="subcounty.blocking")
mod.h7.spill3

###H11: Satisfaction with RS
nod.h11 <- ols(e.satisfaction.b3 ~ treat + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h11

nod.h11.spill <- ols(e.satisfaction.b3 ~ treat + S1 + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h11.spill

nod.h11.spill2 <- ols(e.satisfaction.b3 ~ treat + S2 + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h11.spill2

nod.h11.spill3 <- ols(e.satisfaction.b3 ~ treat + S3 + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h11.spill3

###H10: Satisfaction with Park Management
nod.h10mgmt <- ols(e.satisfaction.b2 ~ treat + b.satisfaction.b2 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h10mgmt

nod.h10mgmt.spill <- ols(e.satisfaction.b2 ~ treat + S1 + b.satisfaction.b2 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h10mgmt.spill

nod.h10mgmt.spill2 <- ols(e.satisfaction.b2 ~ treat + S2 + b.satisfaction.b2 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h10mgmt.spill2

nod.h10mgmt.spill3 <- ols(e.satisfaction.b2 ~ treat + S3 + b.satisfaction.b2 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h10mgmt.spill3

###H13: Support for Conservation (Importance of Protecting Bwindi)
nod.h13 <- ols(e.conservation.b6 ~ treat + b.conservation.b6 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h13

nod.h13.spill <- ols(e.conservation.b6 ~ treat + S1 + b.conservation.b6 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h13.spill

nod.h13.spill2 <- ols(e.conservation.b6 ~ treat + S2 + b.conservation.b6 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h13.spill2

nod.h13.spill3 <- ols(e.conservation.b6 ~ treat + S3 + b.conservation.b6 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking, data=M, cluster="subcounty.blocking")
nod.h13.spill3


######################################################################################
###Table E1. Attrition in endline survey by treatment condition
######################################################################################

table(M$treat,M$attr)
prop.table(table(M$treat,M$attr),1)
chisq.test(table(M$treat,M$attr))


######################################################################################
###Analysis of Attrition reported in-text, SI Section E
######################################################################################
source('outcome_var_setup.R') #Setting up baseline covariates without imputation

#####Analysis on differential missingness reported in text, but not in the table
M.attr <- subset(M, select=c("attr","treat","b.knowledge.b7","b.knowledge.b8","b.knowledge.b9","b.participate.b10","b.participate.b11","b.participate.b12","b.satisfaction.b2","b.satisfaction.b3","b.conservation.b6","subjects.per.village"))
#Note: participation in revenue-sharing meeting is only post-treatment, not included here

#No missingness in any of the main variables
table(M.attr$treat,useNA = "always")
table(M.attr$attr,useNA = "always")
table(M.attr$subjects.per.village,useNA = "always")

M.attr[is.na(M.attr)] <- "missing"

##F-Test on subgroup differences in attrition by baseline survey response
attr.mod <- lm(attr ~ treat + b.knowledge.b7 + b.knowledge.b8 + b.knowledge.b9 + b.participate.b10 + b.participate.b11 + b.participate.b12 + b.satisfaction.b2 + b.satisfaction.b3 + b.conservation.b6 + subjects.per.village, data=M.attr)
summary(attr.mod)

attr.mod.int <- lm(attr ~ treat*b.knowledge.b7 + treat*b.knowledge.b8 + treat*b.knowledge.b9 + treat*b.participate.b10 + treat*b.participate.b11 + treat*b.participate.b12 + treat*b.satisfaction.b2 + treat*b.satisfaction.b3 + treat*b.conservation.b6 + treat*subjects.per.village, data=M.attr)
summary(attr.mod.int)

anova(attr.mod.int,attr.mod) #No indication of differential attrition inconsistent with random chance


######################################################################################
###Table E2. Balance on missingness in pre-treatment measures before and after endline attrition
######################################################################################
M.attr <- subset(M, select=c("attr","treat","b.knowledge.b7","b.knowledge.b8","b.knowledge.b9","b.participate.b10","b.participate.b11","b.participate.b12","b.satisfaction.b2","b.satisfaction.b3","b.conservation.b6","subjects.per.village"))
M.attr[is.na(M.attr)] <- 0
M.attr[M.attr>0] <- 1 #All covariates values transformed to indicator of missinging

###H7.2: Know How RS Works Generally
1 - (sum(M.attr$b.knowledge.b8[M.attr$treat==1]) / length(M.attr$b.knowledge.b8[M.attr$treat==1]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b8[M.attr$treat==1]),
          prob= 1 - (sum(M.attr$b.knowledge.b8[M.attr$treat==1]) / length(M.attr$b.knowledge.b8[M.attr$treat==1]))) / 
     length(M.attr$b.knowledge.b8[M.attr$treat==1])
)

1 - (sum(M.attr$b.knowledge.b8[M.attr$treat==0]) / length(M.attr$b.knowledge.b8[M.attr$treat==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b8[M.attr$treat==0]),
          prob= 1 - (sum(M.attr$b.knowledge.b8[M.attr$treat==0]) / length(M.attr$b.knowledge.b8[M.attr$treat==0]))) / 
     length(M.attr$b.knowledge.b8[M.attr$treat==0])
)

1 - (sum(M.attr$b.knowledge.b8[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.knowledge.b8[M.attr$treat==1 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b8[M.attr$treat==1 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.knowledge.b8[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.knowledge.b8[M.attr$treat==1 & M.attr$attr==0]))) / 
     length(M.attr$b.knowledge.b8[M.attr$treat==1 & M.attr$attr==0])
)

1 - (sum(M.attr$b.knowledge.b8[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.knowledge.b8[M.attr$treat==0 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b8[M.attr$treat==0 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.knowledge.b8[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.knowledge.b8[M.attr$treat==0 & M.attr$attr==0]))) / 
     length(M.attr$b.knowledge.b8[M.attr$treat==0 & M.attr$attr==0])
)

###H7.3: How RS Works in my Village
1 - (sum(M.attr$b.knowledge.b9[M.attr$treat==1]) / length(M.attr$b.knowledge.b9[M.attr$treat==1]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b9[M.attr$treat==1]),
          prob= 1 - (sum(M.attr$b.knowledge.b9[M.attr$treat==1]) / length(M.attr$b.knowledge.b9[M.attr$treat==1]))) / 
     length(M.attr$b.knowledge.b9[M.attr$treat==1])
)

1 - (sum(M.attr$b.knowledge.b9[M.attr$treat==0]) / length(M.attr$b.knowledge.b9[M.attr$treat==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b9[M.attr$treat==0]),
          prob= 1 - (sum(M.attr$b.knowledge.b9[M.attr$treat==0]) / length(M.attr$b.knowledge.b9[M.attr$treat==0]))) / 
     length(M.attr$b.knowledge.b9[M.attr$treat==0])
)

1 - (sum(M.attr$b.knowledge.b9[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.knowledge.b9[M.attr$treat==1 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b9[M.attr$treat==1 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.knowledge.b9[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.knowledge.b9[M.attr$treat==1 & M.attr$attr==0]))) / 
     length(M.attr$b.knowledge.b9[M.attr$treat==1 & M.attr$attr==0])
)

1 - (sum(M.attr$b.knowledge.b9[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.knowledge.b9[M.attr$treat==0 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b9[M.attr$treat==0 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.knowledge.b9[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.knowledge.b9[M.attr$treat==0 & M.attr$attr==0]))) / 
     length(M.attr$b.knowledge.b9[M.attr$treat==0 & M.attr$attr==0])
)

###H9.2: Know Right Person to Contact
1 - (sum(M.attr$b.participate.b12[M.attr$treat==1]) / length(M.attr$b.participate.b12[M.attr$treat==1]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b12[M.attr$treat==1]),
          prob= 1 - (sum(M.attr$b.participate.b12[M.attr$treat==1]) / length(M.attr$b.participate.b12[M.attr$treat==1]))) / 
     length(M.attr$b.participate.b12[M.attr$treat==1])
)

1 - (sum(M.attr$b.participate.b12[M.attr$treat==0]) / length(M.attr$b.participate.b12[M.attr$treat==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b12[M.attr$treat==0]),
          prob= 1 - (sum(M.attr$b.participate.b12[M.attr$treat==0]) / length(M.attr$b.participate.b12[M.attr$treat==0]))) / 
     length(M.attr$b.participate.b12[M.attr$treat==0])
)

1 - (sum(M.attr$b.participate.b12[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.participate.b12[M.attr$treat==1 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b12[M.attr$treat==1 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.participate.b12[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.participate.b12[M.attr$treat==1 & M.attr$attr==0]))) / 
     length(M.attr$b.participate.b12[M.attr$treat==1 & M.attr$attr==0])
)

1 - (sum(M.attr$b.participate.b12[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.participate.b12[M.attr$treat==0 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b12[M.attr$treat==0 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.participate.b12[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.participate.b12[M.attr$treat==0 & M.attr$attr==0]))) / 
     length(M.attr$b.participate.b12[M.attr$treat==0 & M.attr$attr==0])
)

###H9.1: Satisfaction with opportunities to communicate with UWA about RS
1 - (sum(M.attr$b.participate.b11[M.attr$treat==1]) / length(M.attr$b.participate.b11[M.attr$treat==1]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b11[M.attr$treat==1]),
          prob= 1 - (sum(M.attr$b.participate.b11[M.attr$treat==1]) / length(M.attr$b.participate.b11[M.attr$treat==1]))) / 
     length(M.attr$b.participate.b11[M.attr$treat==1])
)

1 - (sum(M.attr$b.participate.b11[M.attr$treat==0]) / length(M.attr$b.participate.b11[M.attr$treat==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b11[M.attr$treat==0]),
          prob= 1 - (sum(M.attr$b.participate.b11[M.attr$treat==0]) / length(M.attr$b.participate.b11[M.attr$treat==0]))) / 
     length(M.attr$b.participate.b11[M.attr$treat==0])
)

1 - (sum(M.attr$b.participate.b11[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.participate.b11[M.attr$treat==1 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b11[M.attr$treat==1 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.participate.b11[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.participate.b11[M.attr$treat==1 & M.attr$attr==0]))) / 
     length(M.attr$b.participate.b11[M.attr$treat==1 & M.attr$attr==0])
)

1 - (sum(M.attr$b.participate.b11[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.participate.b11[M.attr$treat==0 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b11[M.attr$treat==0 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.participate.b11[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.participate.b11[M.attr$treat==0 & M.attr$attr==0]))) / 
     length(M.attr$b.participate.b11[M.attr$treat==0 & M.attr$attr==0])
)

###H8.1: Opportunities to Participate in Planning
1 - (sum(M.attr$b.participate.b10[M.attr$treat==1]) / length(M.attr$b.participate.b10[M.attr$treat==1]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b10[M.attr$treat==1]),
          prob= 1 - (sum(M.attr$b.participate.b10[M.attr$treat==1]) / length(M.attr$b.participate.b10[M.attr$treat==1]))) / 
     length(M.attr$b.participate.b10[M.attr$treat==1])
)

1 - (sum(M.attr$b.participate.b10[M.attr$treat==0]) / length(M.attr$b.participate.b10[M.attr$treat==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b10[M.attr$treat==0]),
          prob= 1 - (sum(M.attr$b.participate.b10[M.attr$treat==0]) / length(M.attr$b.participate.b10[M.attr$treat==0]))) / 
     length(M.attr$b.participate.b10[M.attr$treat==0])
)

1 - (sum(M.attr$b.participate.b10[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.participate.b10[M.attr$treat==1 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b10[M.attr$treat==1 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.participate.b10[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.participate.b10[M.attr$treat==1 & M.attr$attr==0]))) / 
     length(M.attr$b.participate.b10[M.attr$treat==1 & M.attr$attr==0])
)

1 - (sum(M.attr$b.participate.b10[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.participate.b10[M.attr$treat==0 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.participate.b10[M.attr$treat==0 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.participate.b10[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.participate.b10[M.attr$treat==0 & M.attr$attr==0]))) / 
     length(M.attr$b.participate.b10[M.attr$treat==0 & M.attr$attr==0])
)

###H7.1: Satisfaction with RS Info from UWA
1 - (sum(M.attr$b.knowledge.b7[M.attr$treat==1]) / length(M.attr$b.knowledge.b7[M.attr$treat==1]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b7[M.attr$treat==1]),
          prob= 1 - (sum(M.attr$b.knowledge.b7[M.attr$treat==1]) / length(M.attr$b.knowledge.b7[M.attr$treat==1]))) / 
     length(M.attr$b.knowledge.b7[M.attr$treat==1])
)

1 - (sum(M.attr$b.knowledge.b7[M.attr$treat==0]) / length(M.attr$b.knowledge.b7[M.attr$treat==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b7[M.attr$treat==0]),
          prob= 1 - (sum(M.attr$b.knowledge.b7[M.attr$treat==0]) / length(M.attr$b.knowledge.b7[M.attr$treat==0]))) / 
     length(M.attr$b.knowledge.b7[M.attr$treat==0])
)

1 - (sum(M.attr$b.knowledge.b7[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.knowledge.b7[M.attr$treat==1 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b7[M.attr$treat==1 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.knowledge.b7[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.knowledge.b7[M.attr$treat==1 & M.attr$attr==0]))) / 
     length(M.attr$b.knowledge.b7[M.attr$treat==1 & M.attr$attr==0])
)

1 - (sum(M.attr$b.knowledge.b7[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.knowledge.b7[M.attr$treat==0 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.knowledge.b7[M.attr$treat==0 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.knowledge.b7[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.knowledge.b7[M.attr$treat==0 & M.attr$attr==0]))) / 
     length(M.attr$b.knowledge.b7[M.attr$treat==0 & M.attr$attr==0])
)

###H11: Satisfaction with RS
1 - (sum(M.attr$b.satisfaction.b3[M.attr$treat==1]) / length(M.attr$b.satisfaction.b3[M.attr$treat==1]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.satisfaction.b3[M.attr$treat==1]),
          prob= 1 - (sum(M.attr$b.satisfaction.b3[M.attr$treat==1]) / length(M.attr$b.satisfaction.b3[M.attr$treat==1]))) / 
     length(M.attr$b.satisfaction.b3[M.attr$treat==1])
)

1 - (sum(M.attr$b.satisfaction.b3[M.attr$treat==0]) / length(M.attr$b.satisfaction.b3[M.attr$treat==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.satisfaction.b3[M.attr$treat==0]),
          prob= 1 - (sum(M.attr$b.satisfaction.b3[M.attr$treat==0]) / length(M.attr$b.satisfaction.b3[M.attr$treat==0]))) / 
     length(M.attr$b.satisfaction.b3[M.attr$treat==0])
)

1 - (sum(M.attr$b.satisfaction.b3[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.satisfaction.b3[M.attr$treat==1 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.satisfaction.b3[M.attr$treat==1 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.satisfaction.b3[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.satisfaction.b3[M.attr$treat==1 & M.attr$attr==0]))) / 
     length(M.attr$b.satisfaction.b3[M.attr$treat==1 & M.attr$attr==0])
)

1 - (sum(M.attr$b.satisfaction.b3[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.satisfaction.b3[M.attr$treat==0 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.satisfaction.b3[M.attr$treat==0 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.satisfaction.b3[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.satisfaction.b3[M.attr$treat==0 & M.attr$attr==0]))) / 
     length(M.attr$b.satisfaction.b3[M.attr$treat==0 & M.attr$attr==0])
)

###H10: Satisfaction with Park Management
1 - (sum(M.attr$b.satisfaction.b2[M.attr$treat==1]) / length(M.attr$b.satisfaction.b2[M.attr$treat==1]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.satisfaction.b2[M.attr$treat==1]),
          prob= 1 - (sum(M.attr$b.satisfaction.b2[M.attr$treat==1]) / length(M.attr$b.satisfaction.b2[M.attr$treat==1]))) / 
     length(M.attr$b.satisfaction.b2[M.attr$treat==1])
)

1 - (sum(M.attr$b.satisfaction.b2[M.attr$treat==0]) / length(M.attr$b.satisfaction.b2[M.attr$treat==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.satisfaction.b2[M.attr$treat==0]),
          prob= 1 - (sum(M.attr$b.satisfaction.b2[M.attr$treat==0]) / length(M.attr$b.satisfaction.b2[M.attr$treat==0]))) / 
     length(M.attr$b.satisfaction.b2[M.attr$treat==0])
)

1 - (sum(M.attr$b.satisfaction.b2[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.satisfaction.b2[M.attr$treat==1 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.satisfaction.b2[M.attr$treat==1 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.satisfaction.b2[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.satisfaction.b2[M.attr$treat==1 & M.attr$attr==0]))) / 
     length(M.attr$b.satisfaction.b2[M.attr$treat==1 & M.attr$attr==0])
)

1 - (sum(M.attr$b.satisfaction.b2[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.satisfaction.b2[M.attr$treat==0 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.satisfaction.b2[M.attr$treat==0 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.satisfaction.b2[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.satisfaction.b2[M.attr$treat==0 & M.attr$attr==0]))) / 
     length(M.attr$b.satisfaction.b2[M.attr$treat==0 & M.attr$attr==0])
)

###H13: Support for Conservation (Importance of Protecting Bwindi)
1 - (sum(M.attr$b.conservation.b6[M.attr$treat==1]) / length(M.attr$b.conservation.b6[M.attr$treat==1]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.conservation.b6[M.attr$treat==1]),
          prob= 1 - (sum(M.attr$b.conservation.b6[M.attr$treat==1]) / length(M.attr$b.conservation.b6[M.attr$treat==1]))) / 
     length(M.attr$b.conservation.b6[M.attr$treat==1])
)

1 - (sum(M.attr$b.conservation.b6[M.attr$treat==0]) / length(M.attr$b.conservation.b6[M.attr$treat==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.conservation.b6[M.attr$treat==0]),
          prob= 1 - (sum(M.attr$b.conservation.b6[M.attr$treat==0]) / length(M.attr$b.conservation.b6[M.attr$treat==0]))) / 
     length(M.attr$b.conservation.b6[M.attr$treat==0])
)

1 - (sum(M.attr$b.conservation.b6[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.conservation.b6[M.attr$treat==1 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.conservation.b6[M.attr$treat==1 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.conservation.b6[M.attr$treat==1 & M.attr$attr==0]) / length(M.attr$b.conservation.b6[M.attr$treat==1 & M.attr$attr==0]))) / 
     length(M.attr$b.conservation.b6[M.attr$treat==1 & M.attr$attr==0])
)

1 - (sum(M.attr$b.conservation.b6[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.conservation.b6[M.attr$treat==0 & M.attr$attr==0]))
sd(rbinom(n= 10000,
          size= length(M.attr$b.conservation.b6[M.attr$treat==0 & M.attr$attr==0]),
          prob= 1 - (sum(M.attr$b.conservation.b6[M.attr$treat==0 & M.attr$attr==0]) / length(M.attr$b.conservation.b6[M.attr$treat==0 & M.attr$attr==0]))) / 
     length(M.attr$b.conservation.b6[M.attr$treat==0 & M.attr$attr==0])
)


######################################################################################
###Table E3. Balance on values of pre-treatment measures, conditional on non-missingness
######################################################################################
source('outcome_var_setup.R') #Setting up baseline covariates without imputation

M.attr <- subset(M, select=c("attr","treat","b.knowledge.b7","b.knowledge.b8","b.knowledge.b9","b.participate.b10","b.participate.b11","b.participate.b12","b.satisfaction.b2","b.satisfaction.b3","b.conservation.b6","subjects.per.village"))
boot.se <- function(vector,sims=10000){
  store <<- rep(NA, sims)
  for (i in 1:sims){
    store[i] <- mean(sample(vector, replace=T), na.rm=T)
  }
  return(sd(store))
}

###H7.2: Know How RS Works Generally
mean(M.attr$b.knowledge.b8[M.attr$treat==1], na.rm=T)
boot.se(M.attr$b.knowledge.b8[M.attr$treat==1])

mean(M.attr$b.knowledge.b8[M.attr$treat==0], na.rm=T)
boot.se(M.attr$b.knowledge.b8[M.attr$treat==0])

mean(M.attr$b.knowledge.b8[M.attr$treat==1 & M$attr==0], na.rm=T)
boot.se(M.attr$b.knowledge.b8[M.attr$treat==1 & M$attr==0])

mean(M.attr$b.knowledge.b8[M.attr$treat==0 & M$attr==0], na.rm=T)
boot.se(M.attr$b.knowledge.b8[M.attr$treat==0 & M$attr==0])

###H7.3: How RS Works in my Village
mean(M.attr$b.knowledge.b9[M.attr$treat==1], na.rm=T)
boot.se(M.attr$b.knowledge.b9[M.attr$treat==1])

mean(M.attr$b.knowledge.b9[M.attr$treat==0], na.rm=T)
boot.se(M.attr$b.knowledge.b9[M.attr$treat==0])

mean(M.attr$b.knowledge.b9[M.attr$treat==1 & M$attr==0], na.rm=T)
boot.se(M.attr$b.knowledge.b9[M.attr$treat==1 & M$attr==0])

mean(M.attr$b.knowledge.b9[M.attr$treat==0 & M$attr==0], na.rm=T)
boot.se(M.attr$b.knowledge.b9[M.attr$treat==0 & M$attr==0])

###H9.2: Know Right Person to Contact
mean(M.attr$b.participate.b12[M.attr$treat==1], na.rm=T)
boot.se(M.attr$b.participate.b12[M.attr$treat==1])

mean(M.attr$b.participate.b12[M.attr$treat==0], na.rm=T)
boot.se(M.attr$b.participate.b12[M.attr$treat==0])

mean(M.attr$b.participate.b12[M.attr$treat==1 & M$attr==0], na.rm=T)
boot.se(M.attr$b.participate.b12[M.attr$treat==1 & M$attr==0])

mean(M.attr$b.participate.b12[M.attr$treat==0 & M$attr==0], na.rm=T)
boot.se(M.attr$b.participate.b12[M.attr$treat==0 & M$attr==0])

###H9.1: Satisfaction with opportunities to communicate with UWA about RS
mean(M.attr$b.participate.b11[M.attr$treat==1], na.rm=T)
boot.se(M.attr$b.participate.b11[M.attr$treat==1])

mean(M.attr$b.participate.b11[M.attr$treat==0], na.rm=T)
boot.se(M.attr$b.participate.b11[M.attr$treat==0])

mean(M.attr$b.participate.b11[M.attr$treat==1 & M$attr==0], na.rm=T)
boot.se(M.attr$b.participate.b11[M.attr$treat==1 & M$attr==0])

mean(M.attr$b.participate.b11[M.attr$treat==0 & M$attr==0], na.rm=T)
boot.se(M.attr$b.participate.b11[M.attr$treat==0 & M$attr==0])

###H8.1: Opportunities to Participate in Planning
mean(M.attr$b.participate.b10[M.attr$treat==1], na.rm=T)
boot.se(M.attr$b.participate.b10[M.attr$treat==1])

mean(M.attr$b.participate.b10[M.attr$treat==0], na.rm=T)
boot.se(M.attr$b.participate.b10[M.attr$treat==0])

mean(M.attr$b.participate.b10[M.attr$treat==1 & M$attr==0], na.rm=T)
boot.se(M.attr$b.participate.b10[M.attr$treat==1 & M$attr==0])

mean(M.attr$b.participate.b10[M.attr$treat==0 & M$attr==0], na.rm=T)
boot.se(M.attr$b.participate.b10[M.attr$treat==0 & M$attr==0])

###H7.1: Satisfaction with RS Info from UWA
mean(M.attr$b.knowledge.b7[M.attr$treat==1], na.rm=T)
boot.se(M.attr$b.knowledge.b7[M.attr$treat==1])

mean(M.attr$b.knowledge.b7[M.attr$treat==0], na.rm=T)
boot.se(M.attr$b.knowledge.b7[M.attr$treat==0])

mean(M.attr$b.knowledge.b7[M.attr$treat==1 & M$attr==0], na.rm=T)
boot.se(M.attr$b.knowledge.b7[M.attr$treat==1 & M$attr==0])

mean(M.attr$b.knowledge.b7[M.attr$treat==0 & M$attr==0], na.rm=T)
boot.se(M.attr$b.knowledge.b7[M.attr$treat==0 & M$attr==0])

###H11: Satisfaction with RS
mean(M.attr$b.satisfaction.b3[M.attr$treat==1], na.rm=T)
boot.se(M.attr$b.satisfaction.b3[M.attr$treat==1])

mean(M.attr$b.satisfaction.b3[M.attr$treat==0], na.rm=T)
boot.se(M.attr$b.satisfaction.b3[M.attr$treat==0])

mean(M.attr$b.satisfaction.b3[M.attr$treat==1 & M$attr==0], na.rm=T)
boot.se(M.attr$b.satisfaction.b3[M.attr$treat==1 & M$attr==0])

mean(M.attr$b.satisfaction.b3[M.attr$treat==0 & M$attr==0], na.rm=T)
boot.se(M.attr$b.satisfaction.b3[M.attr$treat==0 & M$attr==0])

###H10: Satisfaction with Park Management
mean(M.attr$b.satisfaction.b2[M.attr$treat==1], na.rm=T)
boot.se(M.attr$b.satisfaction.b2[M.attr$treat==1])

mean(M.attr$b.satisfaction.b2[M.attr$treat==0], na.rm=T)
boot.se(M.attr$b.satisfaction.b2[M.attr$treat==0])

mean(M.attr$b.satisfaction.b2[M.attr$treat==1 & M$attr==0], na.rm=T)
boot.se(M.attr$b.satisfaction.b2[M.attr$treat==1 & M$attr==0])

mean(M.attr$b.satisfaction.b2[M.attr$treat==0 & M$attr==0], na.rm=T)
boot.se(M.attr$b.satisfaction.b2[M.attr$treat==0 & M$attr==0])

###H13: Support for Conservation (Importance of Protecting Bwindi)
mean(M.attr$b.conservation.b6[M.attr$treat==1], na.rm=T)
boot.se(M.attr$b.conservation.b6[M.attr$treat==1])

mean(M.attr$b.conservation.b6[M.attr$treat==0], na.rm=T)
boot.se(M.attr$b.conservation.b6[M.attr$treat==0])

mean(M.attr$b.conservation.b6[M.attr$treat==1 & M$attr==0], na.rm=T)
boot.se(M.attr$b.conservation.b6[M.attr$treat==1 & M$attr==0])

mean(M.attr$b.conservation.b6[M.attr$treat==0 & M$attr==0], na.rm=T)
boot.se(M.attr$b.conservation.b6[M.attr$treat==0 & M$attr==0])


######################################################################################
###Figure G1: Reserve power calculation for satisfaction w/ Revenue Sharing
######################################################################################
source('outcome_var_mean_imp_setup.R')
source2('BDD_WD2018_replication.R', start=628, end=629) #sourcing main model results

te.seq <- seq(from=0.01,to=0.5,by=0.01)
power.store <- rep(NA, length(te.seq))

for (i in 1:length(te.seq)){
  nod.h11.power <- felm.ri.power(e.satisfaction.b3 ~ treat + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                                 dta=M.ri, treat.var = "treat", sims=10000, ate.bar = quantile(sort(nod.h11.ri$ate.samp.dist), probs = 0.95), outcome.var = "e.satisfaction.b3", te=te.seq[i])
  power.store[i] <- nod.h11.power$power
}

##Plotting reverse 
par(mar=c(4.1,4.1,1.1,1.1))
plot(x=te.seq, y=power.store, type="l", ylab="Power", xlab="Treatment Effect")
abline(h=0.8, col="red", lty=2)


######################################################################################
###Table G1: Minimum Detectable Positive Effect for Each Survey Item
######################################################################################
source('outcome_var_mean_imp_setup.R')

te.seq <- seq(from=0.2, to=0.8, by=0.01)

###H7.2: Know How RS Works Generally
mod.h7.2.ri <- felm.ri(e.knowledge.b8 ~ treat + b.knowledge.b8 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                       dta=M.ri, treat.var = "treat", sims=10000)

power.store <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  mod.h7.2.power <- lm.ri.power(e.knowledge.b8 ~ treat + b.knowledge.b8 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                                dta=M.ri, treat.var = "treat", sims=10000, ate.bar = quantile(sort(mod.h7.2.ri$ate.samp.dist), probs = 0.95), outcome.var = "e.knowledge.b8", te=te.seq[i])
  power.store[i] <- mod.h7.2.power$power
  if (mod.h7.2.power$power >= 0.8) break
} #ate.bar: 0.08173884
min(te.seq[power.store>=0.8 & !is.na(power.store)]) #0.25

#prop. ceiling obs
table(M$e.knowledge.b8)[length(table(M$e.knowledge.b8))] / sum(table(M$e.knowledge.b8))

#std. effect
0.25/sd(M$e.knowledge.b8[M$treat==0], na.rm=T)

#mde
0.25*(1-0.5111562)


###H7.3: How RS Works in my Village
mod.h7.3.ri <- felm.ri(e.knowledge.b9 ~ treat + b.knowledge.b9 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                       dta=M.ri, treat.var = "treat", sims=10000)

power.store <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  mod.h7.3.power <- lm.ri.power(e.knowledge.b9 ~ treat + b.knowledge.b9 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                                dta=M.ri, treat.var = "treat", sims=10000, ate.bar = quantile(sort(mod.h7.3.ri$ate.samp.dist), probs = 0.95), outcome.var = "e.knowledge.b9", te=te.seq[i])
  power.store[i] <- mod.h7.3.power$power
  if (mod.h7.3.power$power >= 0.8) break
} #ate.bar: 0.07572411
min(te.seq[power.store>=0.8 & !is.na(power.store)]) #0.30

#prop. ceiling obs
table(M$e.knowledge.b9)[length(table(M$e.knowledge.b9))] / sum(table(M$e.knowledge.b9))

#std. effect
0.30/sd(M$e.knowledge.b9[M$treat==0], na.rm=T)

#mde
0.30*(1-0.6135217)


###H9.2: Know Right Person to Contact
mod.h9.2.ri <- lm.ri(e.participate.b12 ~ treat + b.participate.b12 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M.ri, treat.var = "treat", sims=10000)

power.store <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  mod.h9.2.power <- lm.ri.power(e.participate.b12 ~ treat + b.participate.b12 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                                dta=M.ri, treat.var = "treat", sims=10000, ate.bar = quantile(sort(mod.h9.2.ri$ate.samp.dist), probs = 0.95), outcome.var = "e.participate.b12", te=te.seq[i])
  power.store[i] <- mod.h9.2.power$power
  if (mod.h9.2.power$power >= 0.8) break
} #ate.bar: 0.0912536
min(te.seq[power.store>=0.8 & !is.na(power.store)]) #0.6

#prop. ceiling obs
table(M$e.participate.b12)[length(table(M$e.participate.b12))] / sum(table(M$e.participate.b12))

#std. effect
0.60/sd(M$e.participate.b12[M$treat==0], na.rm=T)

#mde
0.60*(1-0.7741607)


###H9.1: Satisfaction with opportunities to communicate with UWA about RS
mod.h9.1.ri <- felm.ri(e.participate.b11 ~ treat + b.participate.b11 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                       dta=M.ri, treat.var = "treat", sims=10000)

power.store <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  mod.h9.1.power <- lm.ri.power(e.participate.b11 ~ treat + b.participate.b11 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                                dta=M.ri, treat.var = "treat", sims=10000, ate.bar = quantile(sort(mod.h9.1.ri$ate.samp.dist), probs = 0.95), outcome.var = "e.participate.b11", te=te.seq[i])
  power.store[i] <- mod.h9.1.power$power
  if (mod.h9.1.power$power >= 0.8) break
} #ate.bar: 0.1017958
min(te.seq[power.store>=0.8 & !is.na(power.store)]) #0.25

#prop. ceiling obs
table(M$e.participate.b11)[length(table(M$e.participate.b11))] / sum(table(M$e.participate.b11))

#std. effect
0.25/sd(M$e.participate.b11[M$treat==0], na.rm=T)

#mde
0.25*(1-0.3679435)


###H8.1: Opportunities to Participate in Planning
mod.h8.1.ri <- felm.ri(e.participate.b10 ~ treat + b.participate.b10 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                       dta=M.ri, treat.var = "treat", sims=10000)

power.store <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  mod.h8.1.power <- lm.ri.power(e.participate.b10 ~ treat + b.participate.b10 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                                dta=M.ri, treat.var = "treat", sims=10000, ate.bar = quantile(sort(mod.h8.1.ri$ate.samp.dist), probs = 0.95), outcome.var = "e.participate.b10", te=te.seq[i])
  power.store[i] <- mod.h8.1.power$power
  if (mod.h8.1.power$power >= 0.8) break
} #ate.bar: 0.1149789
min(te.seq[power.store>=0.8 & !is.na(power.store)]) #0.57

#prop. ceiling obs
table(M$e.participate.b10)[length(table(M$e.participate.b10))] / sum(table(M$e.participate.b10))

#std. effect
0.57/sd(M$e.participate.b10[M$treat==0], na.rm=T)

#mde
0.57*(1-0.7102616)


###H8.3: Participate in RS meetings in past months
M.ri$e.participate.e4 <- as.numeric(M.ri$e.participate.e4)
mod.h8.3.ri <- felm.ri(e.participate.e4 ~ treat + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                       dta=M.ri, treat.var = "treat", sims=10000)

power.store <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  mod.h8.3.power <- lm.ri.power(e.participate.e4 ~ treat + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                                dta=M.ri, treat.var = "treat", sims=10000, ate.bar = quantile(sort(mod.h8.3.ri$ate.samp.dist), probs = 0.95), outcome.var = "e.participate.e4", te=te.seq[i])
  power.store[i] <- mod.h8.3.power$power
  if (mod.h8.3.power$power >= 0.8) break
} #ate.bar: 0.0509953
min(te.seq[power.store>=0.8 & !is.na(power.store)]) #0.29

#prop. ceiling obs
table(M$e.participate.e4)[length(table(M$e.participate.e4))] / sum(table(M$e.participate.e4))

#mde
0.29*(1-0.7377377)


###H7.1: Satisfaction with RS Info from UWA
mod.h7.ri <- felm.ri(e.knowledge.b7 ~ treat + b.knowledge.b7 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                     dta=M.ri, treat.var = "treat", sims=10000)

power.store <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  mod.h7.power <- lm.ri.power(e.knowledge.b7 ~ treat + b.knowledge.b7 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                              dta=M.ri, treat.var = "treat", sims=10000, ate.bar = quantile(sort(mod.h7.ri$ate.samp.dist), probs = 0.95), outcome.var = "e.knowledge.b7", te=te.seq[i])
  power.store[i] <- mod.h7.power$power
  if (mod.h7.power$power >= 0.8) break
} #ate.bar: 0.09747668
min(te.seq[power.store>=0.8 & !is.na(power.store)]) #0.23

#prop. ceiling obs
table(M$e.knowledge.b7)[length(table(M$e.knowledge.b7))] / sum(table(M$e.knowledge.b7))

#std. effect
0.23/sd(M$e.knowledge.b7[M$treat==0], na.rm=T)

#mde
0.23*(1-0.3570712)


###H11: Satisfaction with RS
nod.h11.ri <- felm.ri(e.satisfaction.b3 ~ treat + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                      dta=M.ri, treat.var = "treat", sims=10000)

power.store <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  nod.h11.power <- lm.ri.power(e.satisfaction.b3 ~ treat + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                               dta=M.ri, treat.var = "treat", sims=10000, ate.bar = quantile(sort(nod.h11.ri$ate.samp.dist), probs = 0.95), outcome.var = "e.satisfaction.b3", te=te.seq[i])
  power.store[i] <- nod.h11.power$power
  if (nod.h11.power$power >= 0.8) break
} #ate.bar: 0.1525544
min(te.seq[power.store>=0.8 & !is.na(power.store)]) #0.32

#prop. ceiling obs
table(M$e.satisfaction.b3)[length(table(M$e.satisfaction.b3))] / sum(table(M$e.satisfaction.b3))

#std. effect
0.32/sd(M$e.satisfaction.b2[M$treat==0], na.rm=T)

#mde
0.32*(1-0.2862903)


###H10: Satisfaction with Park Management
nod.h10mgmt.ri <- felm.ri(e.satisfaction.b2 ~ treat + b.satisfaction.b2 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                          dta=M.ri, treat.var = "treat", sims=10000)

power.store <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  nod.h10mgmt.power <- lm.ri.power(e.satisfaction.b2 ~ treat + b.satisfaction.b2 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                                   dta=M.ri, treat.var = "treat", sims=10000, ate.bar = quantile(sort(nod.h10mgmt.ri$ate.samp.dist), probs = 0.95), outcome.var = "e.satisfaction.b2", te=te.seq[i])
  power.store[i] <- nod.h10mgmt.power$power
  if (nod.h10mgmt.power$power >= 0.8) break
} #ate.bar: 0.1390786
min(te.seq[power.store>=0.8 & !is.na(power.store)]) #0.37

#prop. ceiling obs
table(M$e.satisfaction.b2)[length(table(M$e.satisfaction.b2))] / sum(table(M$e.satisfaction.b2))

#std. effect
0.37/sd(M$e.satisfaction.b2[M$treat==0], na.rm=T)

#mde
0.37*(1-0.4407631)


###H13: Support for Conservation (Importance of Protecting Bwindi)
nod.h13.ri <- lm.ri(e.conservation.b6 ~ treat + b.conservation.b6 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                    dta=M.ri, treat.var = "treat", sims=10000)

power.store <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  nod.h13.power <- lm.ri.power(e.conservation.b6 ~ treat + b.conservation.b6 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                               dta=M.ri, treat.var = "treat", sims=10000, ate.bar = quantile(sort(nod.h13.ri$ate.samp.dist), probs = 0.95), outcome.var = "e.conservation.b6", te=te.seq[i])
  power.store[i] <- nod.h13.power$power
  if (nod.h13.power$power >= 0.8) break
} #ate.bar: 0.03871926
min(te.seq[power.store>=0.8 & !is.na(power.store)]) #0.60

#prop. ceiling obs
table(M$e.conservation.b6)[length(table(M$e.conservation.b6))] / sum(table(M$e.conservation.b6))

#std. effect
0.60/sd(M$e.conservation.b6[M$treat==0], na.rm=T)

#mde
0.60*(1-0.9078156)

#for paragraph
table(M$e.conservation.b6)


######################################################################################
###Table G2: Treatment Effects on Indices of Knowledge, Participation, and Satisfaction
######################################################################################
source('outcome_var_mean_imp_setup.R')

M$e.participate.e4 <- as.numeric(M$e.participate.e4)

M$e.knowledge <- M$e.knowledge.b8*1.25 + M$e.knowledge.b9*1.25 + M$e.participate.b12*1.25
M$e.participate <- M$e.participate.b10*1.25 + M$e.participate.e4*5 + M$e.participate.b11
M$e.satisfaction <- M$e.knowledge.b7 + M$e.satisfaction.b2 + M$e.satisfaction.b3 + M$e.conservation.b6*1.25
#Note: NA values on single items will result in NA values on index automatically

#To find out proportion at ceiling
table(M$e.knowledge)
table(M$e.participate)
table(M$e.satisfaction)

M$b.knowledge <- M$b.knowledge.b8*1.25 + M$b.knowledge.b9*1.25 + M$b.participate.b12*1.25
M$b.participate <- M$b.participate.b10*1.25 + M$b.participate.b11
M$b.satisfaction <- M$b.knowledge.b7 + M$b.satisfaction.b2 + M$b.satisfaction.b3 + M$b.conservation.b6*1.25

M.ri <- merge(M, sbj_RI, by="Subject_ID_number", all.x=TRUE)

###RI inference on the indices, generating ate.bar from felm.ri

knowledge.ri <- felm.ri(e.knowledge ~ treat + b.knowledge + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                        dta=M.ri, treat.var = "treat", sims=10000)
#ate: -0.005619
#se: 0.1484776
#p.two.way: 0.9709
#p.one.way.greater: 0.5247
#p.one.way.lesser: 0.4753
#quantile(sort(knowledge.ri$ate.samp.dist), probs = 0.95) #0.2481223

participate.ri <- felm.ri(e.participate ~ treat + b.participate + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                          dta=M.ri, treat.var = "treat", sims=10000)
#ate: -0.3809
#se: 0.2520316
#p.two.way: 0.1325
#p.one.way.greater: 0.9356
#p.one.way.lesser: 0.0644
#quantile(sort(participate.ri$ate.samp.dist), probs = 0.95) #0.415441

satisfaction.ri <- felm.ri(e.satisfaction ~ treat + b.satisfaction + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                           dta=M.ri, treat.var = "treat", sims=10000)
#ate: 0.012490
#se: 0.2210442
#p.two.way: 0.9544
#p.one.way.greater: 0.4812
#p.one.way.lesser: 0.5188
#quantile(sort(satisfaction.ri$ate.samp.dist), probs = 0.95) #0.3641675


######################################################################################
###Figure G2: Power on Family-Based Indices
#Note: must run code for Table G2
######################################################################################
#Note: lm.ri.power is faster (no error clustering), but gives same results

te.seq <- seq(from=0.05,to=0.95,by=0.05)

power.store.k <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  k.power <- lm.ri.power(e.knowledge ~ treat + b.knowledge + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                         dta=M.ri, treat.var = "treat", sims=10000, ate.bar = 0.2481223, outcome.var = "e.knowledge", te=te.seq[i])
  power.store.k[i] <- k.power$power
}

power.store.p <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  p.power <- lm.ri.power(e.participate ~ treat + b.participate + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                         dta=M.ri, treat.var = "treat", sims=10000, ate.bar = 0.415441, outcome.var = "e.participate", te=te.seq[i])
  power.store.p[i] <- p.power$power
}

power.store.s <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  s.power <- lm.ri.power(e.satisfaction ~ treat + b.satisfaction + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                         dta=M.ri, treat.var = "treat", sims=10000, ate.bar = 0.3641675, outcome.var = "e.satisfaction", te=te.seq[i])
  power.store.s[i] <- s.power$power
}

min(te.seq[power.store.k>=0.8]) #0.65
min(te.seq[power.store.p>=0.8]) #0.85
min(te.seq[power.store.s>=0.8]) #0.65

##Plotting reverse 

#pdf(file="FigH2_IndexPower_170907.pdf", width=6.5, height=5)
par(mar=c(3.6,3.1,1.1,1.1))
par(mgp=c(2,1,0))
par(oma=c(0,0,0,0))
par(mfrow=c(3,1))

plot(x=te.seq, y=power.store.k, type="l", ylab="Power", xlab="Treatment Effect", main="Knowledge")
abline(h=0.8, col="red", lty=2)

plot(x=te.seq, y=power.store.p, type="l", ylab="Power", xlab="Treatment Effect", main="Participation")
abline(h=0.8, col="red", lty=2)

plot(x=te.seq, y=power.store.s, type="l", ylab="Power", xlab="Treatment Effect", main="Satisfaction")
abline(h=0.8, col="red", lty=2)
#dev.off()

##Precise power calculations for MDE

te.seq <- seq(from=0.6,to=0.65,by=0.001)
power.store.k <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  k.power <- lm.ri.power(e.knowledge ~ treat + b.knowledge + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                         dta=M.ri, treat.var = "treat", sims=10000, ate.bar = 0.2481223, outcome.var = "e.knowledge", te=te.seq[i])
  power.store.k[i] <- k.power$power
}

te.seq <- seq(from=0.8,to=0.85,by=0.001)
power.store.p <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  p.power <- lm.ri.power(e.participate ~ treat + b.participate + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                         dta=M.ri, treat.var = "treat", sims=10000, ate.bar = 0.415441, outcome.var = "e.participate", te=te.seq[i])
  power.store.p[i] <- p.power$power
}

te.seq <- seq(from=0.6,to=0.65,by=0.001)
power.store.s <- rep(NA, length(te.seq))
for (i in 1:length(te.seq)){
  s.power <- lm.ri.power(e.satisfaction ~ treat + b.satisfaction + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                         dta=M.ri, treat.var = "treat", sims=10000, ate.bar = 0.3641675, outcome.var = "e.satisfaction", te=te.seq[i])
  power.store.s[i] <- s.power$power
}

min(te.seq[power.store.k>=0.8]) #0.603
min(te.seq[power.store.s>=0.8]) #0.631

te.seq <- seq(from=0.8,to=0.85,by=0.001)
min(te.seq[power.store.p>=0.8]) #0.845

mean(M$e.knowledge, na.rm=T)


######################################################################################
###Table H1: Modeling number of subjects per village as a function of village-level characteristics
######################################################################################
source('outcome_var_setup.R')

pull.mean <- function(source, var, id){
  id.list <- sort(unique(source[,id]))
  store <- rep(NA, length(unique(source[,id])))
  for (i in 1:length(unique(source[,id]))){
    store[i] <- mean(source[source[,id]==id.list[i],var], na.rm=TRUE)
  }
  return(store)
}

village.dta <- data.frame(sort(unique(M$Updated.Village)))
village.dta$subjects.per.village <- pull.mean(source=M, var="subjects.per.village", id="Updated.Village")
village.dta$b.knowledge <- pull.mean(source=M, var="b.knowledge", id="Updated.Village")
village.dta$b.participate <- pull.mean(source=M, var="b.participate", id="Updated.Village")
village.dta$b.satisfaction <- pull.mean(source=M, var="b.satisfaction", id="Updated.Village")
village.dta$age <- pull.mean(source=M, var="E_Subject_s_age", id="Updated.Village")

M$female <- ifelse(M$E_subject_s_gender=="female",1,0)
village.dta$female <- pull.mean(source=M, var="female", id="Updated.Village")

M$literate <- ifelse(M$E_Can_you_read_001=="yes__i_can_rea",1,0) #total literacy
village.dta$literate <- pull.mean(source=M, var="literate", id="Updated.Village")

M$poverty <- ifelse(M$E_Subject_s_income=="20_000_shillin" | M$E_Subject_s_income=="20_000___100_0", 1, 0)
M$poverty[M$E_Subject_s_income=="refused_to_ans"] <- NA
village.dta$poverty <- pull.mean(source=M, var="poverty", id="Updated.Village")

mod <- lm(subjects.per.village ~ b.knowledge + b.participate + b.satisfaction + age + female + literate + poverty, data=village.dta)
summary(mod)


######################################################################################
###Table H2: Modeling number of subjects per village as a function of village-level characteristics
######################################################################################
source('outcome_var_mean_imp_setup.R')

###H7.2: Know How RS Works Generally
mod.h7.2.d <- lm(e.knowledge.b8 ~ treat*subjects.per.village + b.knowledge.b8 + subcounty.blocking, data=M)
summary(mod.h7.2.d)

###H7.3: How RS Works in my Village
mod.h7.3.d <- lm(e.knowledge.b9 ~ treat*subjects.per.village + b.knowledge.b9 + subcounty.blocking, data=M)
summary(mod.h7.3.d)

###H9.2: Know Right Person to Contact
mod.h9.2.d <- lm(e.participate.b12 ~ treat*subjects.per.village + b.participate.b12 + subcounty.blocking, data=M)
summary(mod.h9.2.d)

###H9.1: Satisfaction with opportunities to communicate with UWA about RS
mod.h9.1.d <- lm(e.participate.b11 ~ treat*subjects.per.village + b.participate.b11 + subcounty.blocking, data=M)
summary(mod.h9.1.d)

###H8.1: Opportunities to Participate in Planning
mod.h8.1.d <- lm(e.participate.b10 ~ treat*subjects.per.village + b.participate.b10 + subcounty.blocking, data=M)
summary(mod.h8.1.d)

###H8.3: Participate in RS meetings in past months
mod.h8.3.d <- lm(e.participate.e4 ~ treat*subjects.per.village + subcounty.blocking, data=M)
summary(mod.h8.3.d)

###H7.1: Satisfaction with RS Info from UWA
mod.h7.d <- lm(e.knowledge.b7 ~ treat*subjects.per.village + b.knowledge.b7 + subcounty.blocking, data=M)
summary(mod.h7.d)

###H11: Satisfaction with RS
nod.h11.d <- lm(e.satisfaction.b3 ~ treat*subjects.per.village + b.satisfaction.b3 + subcounty.blocking, data=M)
summary(nod.h11.d)

###H10: Satisfaction with Park Management
nod.h10mgmt.d <- lm(e.satisfaction.b2 ~ treat*subjects.per.village + b.satisfaction.b2 + subcounty.blocking, data=M)
summary(nod.h10mgmt.d)

###H13: Support for Conservation (Importance of Protecting Bwindi)
nod.h13.d <- lm(e.conservation.b6 ~ treat*subjects.per.village + b.conservation.b6 + subcounty.blocking, data=M)
summary(nod.h13.d)

